── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.2 ✔ tibble 3.3.0
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.1.0
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Code
library(sf)
Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
Code
library(tmap)library(terra)
terra 1.8.60
Attaching package: 'terra'
The following object is masked from 'package:tidyr':
extract
Code
library(stars)
Loading required package: abind
Read data
Code
# Read night lights datatile05_07 <-read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021038.h08v05.001.2021039064328.tif'))tile06_07 <-read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021038.h08v06.001.2021039064329.tif'))tile05_16 <-read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather','data', 'VNP46A1', 'VNP46A1.A2021047.h08v05.001.2021048091106.tif'))tile06_16 <-read_stars(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'VNP46A1', 'VNP46A1.A2021047.h08v06.001.2021048091105.tif'))# Read roads datahighways <-st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'gis_osm_roads_free_1.gpkg'),query ="SELECT * FROM gis_osm_roads_free_1 WHERE fclass = 'motorway'") %>%st_transform(crs ='EPSG:3083')# Read houses datahouses <-st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data','gis_osm_buildings_a_free_1.gpkg'),query ="SELECT * FROM gis_osm_buildings_a_free_1 WHERE (type IS NULL AND name IS NULL) OR type IN ('residential', 'apartments', 'house', 'static_caravan', 'detached')") %>%st_transform('EPSG:3083')# Explore the contents of the geodatabasesocioeconomic <-st_layers(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'))# Extract the census tract informationcensus_tract <-st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'), layer ='ACS_2019_5YR_TRACT_48_TEXAS') %>%st_transform('EPSG:3083')# Extract the income informationincome <-st_read(here::here('posts', '2025-12-impacts-of-extreme-weather', 'data', 'ACS_2019_5YR_TRACT_48_TEXAS.gdb'), layer ='X19_INCOME')
Merge raster data and create houston bounding box
Code
# Merge two raster data of February 7nightlight_07 <-st_mosaic(tile05_07, tile06_07)# Merge two raster data of February 16nightlight_16 <-st_mosaic(tile05_16, tile06_16)# Create Houston bounding boxhouston_bbox <-st_bbox(c(xmin =-96.5, xmax =-94.5, ymin =29, ymax =30.5),crs ='EPSG:4326')# Crop all raster data to Houston areanightlight_07_crop <-st_crop(nightlight_07, houston_bbox)nightlight_07_crop[(nightlight_07_crop >1000) | (nightlight_07_crop <=0)] <-NAnightlight_16_crop <-st_crop(nightlight_16, houston_bbox)nightlight_16_crop[(nightlight_16_crop >1000) | (nightlight_16_crop <=0)] <-NA
Visualize before and after night light intensities
# Calculate the difference in light intensitynightlight_diff <- nightlight_07 - nightlight_16# Crop and reclassify the raster data mask <-st_crop(nightlight_diff, houston_bbox)mask[mask <200] <-NA# Vectorize the blackout maskblackout <-st_as_sf(mask,as_points =FALSE,merge =TRUE) %>%st_make_valid() %>%st_transform(crs ='EPSG:3083')
Code
# Combine all highway geometries into onehighways_union <-st_union(highways)# Create 200m buffer around all highwayshighways_buffer <-st_buffer(highways_union, dist =200)# Find areas that experienced blackouts further than 200m from highwaysblackout_far <-st_difference(blackout, highways_buffer)
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Code
# Visualize the blackout in Houstontm_basemap('OpenStreetMap') +tm_shape(blackout_far) +tm_polygons(fill ='ivory', fill_alpha =0.8, col ='black') +tm_title('Blackout areas in Houston', fontface ='bold', size =1) +tm_layout(legend.show =TRUE) +tm_scalebar(position =c('left', 'bottom')) +tm_compass(position =c('right', 'top'))
Estimate of the number of homes in Houston that lost power
Approximately 157,970 homes in Houston lost power in February 2021 due to the severe winter storms and resulting major power crisis.
Note: I used st_intersects. st_intersects returns a boolean value indicating whether two geometries share any features or if they have any point in common. TRUE (or 1) if the geometries intersect, FALSE (or 0) if they are completely separate.
Code
# Keep buildings that intersect with blackout areashouses_blackout <-st_intersects(houses, blackout_far)# Check if each building is in a blackout areain_blackout <-lengths(houses_blackout) >0# Filter buildings in blackout areasblackout_houses <- houses[in_blackout, ]# Number of impacted buildingsn_blackout_houses <-nrow(blackout_houses)# Convert all buildings to centroidshouses_points <- houses %>%st_centroid()
Warning: st_centroid assumes attributes are constant over geometries
Code
# Houses impacted by blackoutshouses_impacted <- blackout_houses %>%st_centroid()
Warning: st_centroid assumes attributes are constant over geometries
Code
# Houses not impacted by blackoutshouses_not_impacted <- houses_points %>%filter(!osm_id %in% blackout_houses$osm_id)
# Rename the column in income dataset to match with the census tractincome_renamed <- income %>%rename("GEOID_Data"="GEOID")# Left join the two datasets by the geometrycensus_income <- census_tract %>%left_join(income_renamed, by ='GEOID_Data')
Code
# Check CRS of the datasetif(st_crs(census_income) ==st_crs(houston_bbox)){print("Coordinate reference systems match!")} else {warning("Update coordinate reference systems to match!")}
Warning: Update coordinate reference systems to match!
Code
# Transform the CRS and crop to houston bounding boxcensus_income <-st_transform(census_income, crs ='epsg:4326') # Crop census income to Houstoncensus_income_crop <- census_income %>%st_crop(houston_bbox)
Warning: attribute variables are assumed to be spatially constant throughout
all geometries
Code
# Transform and clean blackout housesblackout_houses_transformed <- blackout_houses %>%st_transform(crs ='epsg:4326')# Identify the blackout areascensus_blackout <-st_intersects(census_income_crop, blackout_houses_transformed)# Create the columncensus_blackout_col <- census_income_crop %>%mutate(blackout =lengths(census_blackout) >0)# Select only the values that are true for plottingcensus_blackout_col_filtered <- census_blackout_col[census_blackout_col$blackout, ]
Map of the census tracts in Houston that lost power
Code
# Map the census tracts that lost powertm_basemap('CartoDB.VoyagerNoLabels') +tm_shape(census_blackout_col_filtered) +tm_polygons(fill ='blackout',fill.scale =tm_scale_categorical(values =c('TRUE'='firebrick'),labels =c('TRUE'='Blackout')),fill.legend =tm_legend(title ='Blackout status'),col ='grey50',lwd =0.2) +tm_title('Census tracts in Houston that lost power',position =tm_pos_out('center', 'top'),fontface ='bold',size =1) +tm_layout(legend.outside =TRUE,legend.outside.position ='right',legend.bg.color ='white',legend.bg.alpha =0.8) +tm_compass(position =c('right', 'bottom'),type ='arrow',size =1) +tm_scalebar(position =c('left', 'top'),text.size =0.8)
Plot comparing the distributions of median household income
Code
ggplot(census_blackout_col %>%filter(!is.na(B19013e1)), # Remove NA valuesaes( x = blackout, y = B19013e1, fill = blackout)) +geom_boxplot(alpha =0.8) +scale_x_discrete(labels =c('TRUE'='Blackout','FALSE'='No blackout')) +scale_fill_manual(values =c('TRUE'='firebrick','FALSE'='midnightblue'),labels =c('Blackout', 'No blackout')) +labs(title ='Median household income for census tracts',subtitle ='Blackout status in Houston census tracts',x ='Blackout status',y ='Median household income in $',fill ='Status') +scale_y_continuous(labels = scales::dollar_format()) +theme_minimal() +theme(plot.title =element_text(face ='bold', size =16),plot.subtitle =element_text(size =14),legend.position ='bottom')
Figure 1: This plot compares the distributions of median household income between census tracts that experienced blackouts and those that did not during the severe winter storms of February 2021 in Houston, Texas.
Brief reflection
A brief reflection (approx. 100 words) summarizing your results and discussing any limitations to this study
Answer This analysis used satellite data to compare median household incomes between areas affected and unaffected by the blackout. Interestingly, higher-income census tracts were more affected, as their median income was slightly higher than that of non-blackout areas. This result may reflect several limitations such as wealthier neighborhoods often have more outdoor lighting, census tract–level data may overlook neighborhood differences, and blackout areas could be misclassified due to highway lighting or cloud cover. The analysis may also have assumed all structures were residential, while vegetation near power lines and the 200-meter buffer used could further limit accuracy of true lighting.
Citation
BibTeX citation:
@online{poudel2025,
author = {Poudel, Aakriti},
title = {Identifying the Impacts of Extreme Weather},
date = {2025-12-03},
url = {https://aakriti-poudel-chhetri.github.io/posts/2025-12-impacts-of-extreme-weather/},
langid = {en}
}